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In this study, we investigate the phenomenon of collective motion in binary mixtures of self- 
propelled particles. More precisely, we consider two particle species, each of which consisting of 
pointlike objects that propel with a velocity of constant magnitude. Within each species, the 
particles try to achieve polar alignment of their velocity vectors, whereas we analyze the cases of 
preferred polar, antiparallel, as well as perpendicular alignment between particles of different species. 
Our focus is on the effect that the interplay between the two species has on the threshold densities 
for the onset of collective motion and on the nature of the solutions above onset. For this purpose, 
we start from suitable Langevin equations in the particle picture, from which we derive mean field 
equations of the Fokker-Planck type and finally macroscopic continuum field equations. We perform 
particle simulations of the Langevin equations, linear stability analyses of the Fokker-Planck and 
macroscopic continuum equations, and we numerically solve the Fokker-Planck equations. Both, 
spatially homogeneous and inhomogeneous solutions are investigated, where the latter correspond 
to stripe-like flocks of collectively moving particles. In general, the interaction between the two 
species reduces the threshold density for the onset of collective motion of each species. However, 
this interaction also reduces the spatial organization in the stripe-like flocks. The case that shows the 
most interesting behavior is the one of preferred perpendicular alignment between different species. 
There, a competition between polar and truly nematic orient at ional ordering of the velocity vectors 
takes place within each particle species. Finally, depending on the alignment rule for particles of 
different species and within certain ranges of particle densities, identical and inverted spatial density 
profiles can be found for the two particle species. The system under investigation is confined to two 
spatial dimensions. 

PACS numbers: 87.18. Gh, 64.75. Gh, 05.10.Gg, 64.75.Cd 



I. INTRODUCTION 

Inspired by the biological observation of microorgan- 
isms, self-propulsion has widely been studied as the mo- 
tion of microswimmers through their viscous fluid envi- 
ronment. Different models were suggested and analyzed 
for this kind of self-propulsion in the low Reynolds num- 
ber limit. Examples are active rods propagating planar 
or spiral waves along their body [1] |2], Purcell's swim- 
mer that consists of three bars connected by two joints 
[H H] , a model swimmer of three linked spheres [5 , or 
helical filaments that propagate kinks between regions of 
opposite handedness through their body 

Apart from the propulsion mechanism of single iso- 
lated objects, their collective behavior under hydrody- 
namic coupling has been moved into the focus of recent 
investigations. For example, hydrodynamic continuum 
equations have been derived by coarse-graining the parti- 
cle picture of interacting microswimmers [7 . In this case, 
orientational order parameters follow from locally aver- 
aging the single particle velocities. The mode structure 
and hydrodynamic instabilities of the corresponding con- 
tinuum equations were analyzed [THS]. Here, we can iden- 
tify the thresholds at which orientationally non-ordered 
states become unstable with the onset of collective mo- 
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tion. Aspects of synchronization of active particle mo- 
tion mediated by hydrodynamic interactions have been 
studied in detail p^HIS] , Furthermore, a recent work 
incorporates the effect of an axial macroscopic dynamic 
variable into the hydrodynamic framework [M]. Such 
a macroscopic non-equilibrium variable may appear, for 
example, in the case of collectively rotating helical fila- 
ments. 

In this paper, we will restrict ourselves to collective 
two-dimensional motion of self-propelled particles on a 
substrate. Therefore, no long-range hydrodynamic in- 
teractions are taken into account. Such a system was 
used, for example, to study the motion of deformable 
self-propelled particles and domains p!5Hl8] . There, the 
velocity of each particle was determined through cou- 
pling to its deformation (and vice versa). In contrary, 
in simpler models, a constant driving force is directly 
added to the momentum equation for each particle [19]- 
(2T]. It balances the friction forces acting on the particle. 
The driving force is oriented, for example, along the long 
axis of rod-like particles [20l Hi] . Shape dependent short- 
range interactions as, e.g., excluded volume interactions 
can lead to alignment and thus collective motion. Again, 
macroscopic equations were derived and investigated for 
these cases [21 . 

When we are only interested in the basic mechanism 
that leads to collective motion, we only keep the basic 
necessary features in a minimal model. Then the velocity 
magnitude for each particle can simply be kept fixed as 
a constant in time [22 . Only the position and the angle 
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characterizing its direction of motion are maintained as 
variables for each particle. For such a system, we can 
define local alignment rules for the velocity vectors of two 
interacting particles [23-28 . Since momentum exchange 
with the substrate is possible, the total momentum of the 
interacting particles does not need to be conserved. 

Such approaches can be seen as variants of the model 
introduced by Vicsek et al. [29l [30] and were extensively 
studied numerically. Both, polar [23l [25H271 [3TH33] and 
nematic [24 , 28 , 32 , 33 alignment rules were investigated. 
In the first case, velocity vectors in the ordered state 
point into the same direction, whereas in the second case 
antiparallel alignment of velocity vectors is equally al- 
lowed. The focus was on the transition from disordered 
to collective, i.e. macroscopically ordered, motion. This 
transition occurs with decreasing orientational noise or 
increasing particle density and was found to be discontin- 
uous for sufficiently large system sizes [23l [31] . An inter- 
esting issue was the emergence of spatial heterogeneities 
above (but close to) threshold [23 } 124 } [28 } [3Qti33] . More 
precisely, ffocks and stripes of high particle density and 
orientational order were observed. They could move with 
velocity magnitudes close to the single particle speed. 

Apart from that, hydrodynamic continuum equations 
were derived from alignment rules through a Boltz- 
mann approach [25l[26]. Furthermore, continuum equa- 
tions of the Fokker-Planck type were obtained using the 
mean field approximation [34l [35j. The features listed 
above could be reproduced in numerical investigations of 
macroscopic continuum equations [36l [37] . 

We will follow in this paper the outlined minimal model 
approach. The system that we will study is a binary mix- 
ture of self-propelled particles. Previously, systems of dif- 
ferent self-propelled particles were analyzed in predator- 
prey scenarios [38l [39]. In the current work, however, 
we focus on the effect of alignment interactions between 
different particle species. More precisely, we concentrate 
on the onset and features of collective motion that re- 
sult from the interaction between two different groups of 
particles. Therefore, in our analysis, the main parameter 
of interest will be the orientational coupling parameter 
between particles of different species. 

In the next section, we introduce in detail our minimal 
model of binary self-propelled particle mixtures on the 
particle level. From these, we derive mean field equa- 
tions in the spirit of the Fokker-Planck approach. After 



that, in section III, these equations are checked for lin- 
ear instabilities of orientational order as a function of 
the averaged particle densities. We study the cases of 
polar and antiparallel alignment of the velocity vectors 
as well as a preferred perpendicular alignment. These 
cases were further investigated by particle simulations, 
the results of which are presented in section IV After 



that, we compare to numerical results obtained from the 
Fokker-Planck approach in section [Vj Finally, we derive 
macroscopic continuum equations in section [VI[ and dis- 
cuss the stability of their solutions, before we conclude. 



II. PARTICLE AND FIELD DESCRIPTION 

In the following, we describe our two-dimensional min- 
imal model of a binary mixture of self-propelled particles. 
We consider N interacting particles. From these N par- 
ticles, Ni particles are of species 1, the other N — Ni 
particles are of species 2. Each particle is assumed to 
propel with a velocity of fixed constant magnitude. More 
precisely, particles of species 1 and 2 propel with a ve- 
locity of constant magnitude 2ui and 2u2^ respectively. 
The orientation of the velocity vector of each particle i 
(i = 1, 2, . . . , A^) can be characterized by a single orienta- 
tion angle Oi in the two-dimensional plane. We measure 
these angles Oi with respect to the positive x-axis. The 
independent variables in our model are thus the position 
vectors and the velocity orientation angles Oi of all 
particles, i = 1, 2, . . . , A/". 

These variables follow a set of coupled Langevin equa- 
tions 



dvi 
dt 

dO^ 
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Vi{0i) = 2u 
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dU{ri,...,rN,Oi,...,ON) 
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Here, i^(^) = Ui for i = 1,...,^"! and i^(^) = U2 for 
i = A"! + 1, . . . , A". ri(t) is a Gaussian stochastic force of 
zero mean and correlation of the form (r^(t)rj(t')) = 
2DiSij5{t — t'). Likewise, we have Di = Di for i = 
1, . . . , A^i and A = ^2 for z = A"! + 1, . . . , A". 

In this context, /7(ri, . . . , r^r, 6>i, • • • , 6>Ar) plays the role 
of an interaction potential between the particles. It con- 
trols the preferred orientational alignment (or misalign- 
ment) between particles of identical and different species. 
We assume the following form of a pairwise particle- 
particle interaction potential, 

[/(ri, . . . ,rAr,6>i, . . .,0n) 
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(3) 



Here, we assume that > and ^2 > 0. The terms with 
the coefficients gi and ^2 make two particles of the same 
species try to align their velocity vectors parallel to each 
other and pointing into the same direction (polar align- 
ment). In this study, we restrict ourselves to idealized 
local point-like interactions, which leads to the spatial 
5- functions in Eq. (|3|. 

For particles of different species, we consider the three 
different alignment rules depicted in Fig. [l] Parallel 
alignment of the velocity vectors is achieved through the 
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antiparallel 



perpendicular 



FIG. 1. Illustration of polar, antiparallel, and perpendicu- 
lar alignment interactions between particles of the different 
species 1 and 2. 



angular interaction potential 

Ug{0, - Oj) = Ulie, - Oj) = cos{0, - Oj). (4) 

We find preferred polar alignment for ^ > and preferred 
antiparallel alignment for g < 0. Due to the idealized 
local point-like interactions, antiparallel alignment does 
not lead to spatial repulsion. Flocks of different particle 
species can simply penetrate each other. In this case, for 
the two species to spatially avoid each other, spatially 
nonlocal repulsive forces and/or excluded volume forces 
must be taken into account. 

In the idealized picture, for particles of different species 
to try to spatially avoid each other, the preferred align- 
ment of their velocity vectors must not be parallel to each 
other. We investigate preferred perpendicular alignment, 
which is supported by the term with the coefficient g via 

UgiOi- ) = {0i -0j) = sin^ {0i - Oj ) (5) 

and ^ > 0. 

We now ask for a temporal evolution equation for the 
probability density /(ri, . . . , r^r, ^i, . . . , ^at, of finding 
simultaneously particles 1, . . . , A/" at positions ri, . . . , tat 
and with velocity orientations . . . ,0^- Through the 
well-known procedures [40l [41] we obtain an equation of 
the Fokker-Planck type that reads 

9/(ri,...,rAr,6>i,...,6>Ar,t) 
dt 

N . 

= ^ ^ — "^Ui [ cos Oi dx + sin Oi dy\ 

S[/(ri, . . . ,rAr,6>i, . . . ,6>Ar) 

^^'^ m 

+ Aa,',|/(ri,...,r^,^i,...,^^,t). (6) 



From this equation, we find the time evolution of the 
one-particle density p^"^^(r, 6>, by integrating out all 
variables r2, . . . , r^r, ^2, • • • , ^at, identifying the remain- 
ing variables ri = r and 9i = 0^ and multiplying by Ni. 
Here, the superscript of p^^^ denotes the one-particle den- 
sity, whereas the subscript refers to species 1. Likewise, 
we derive a dynamic equation for the one-particle density 
P2"^^(r, of species 2. 

These equations contain the two-particle densities. 

(2) 

For example, the two-particle density p\-{ {r^r^O^O^) 
is obtained from /(ri, . . . , tat, 6>i, . . . , ^at, ^) by inte- 
grating out all variables rs , . . . , r at , ^3 , . . . , ^at , iden- 
tifying the remaining variables ri = r, r2 = r', 
61 = 0^ as well as 62 = and multiplying by 
A^i(A^i — 1). To close the equations, we apply the 
mean field approximation for the two-particle densities 

p['^ (r, 0)pi'^ (r^ ^0, P^i^ (r, r^ 0, 0') = pi'^ (r, 0)p['^ (r^ ^0 , 
and p^2^(r,r^6>,6>0 = pi^\r,0)p'^^^\r' ,0'). The analogous 
procedure was applied before in the case of a single parti- 
cle species [25 , 26 , 35 . In ref. [42 , a similar procedure is 
applied to the conventional case of non-propelled Brow- 
nian particle mixtures. 

Since we are only referring to one-particle densities, 
we omit the superscript during the rest of this 
manuscript. The resulting equations read 

dpi{r,0,t) ^ 
dt 

— 2ui [ cos dx -\- sin 6 dy] pi (r, 0, t) 

^Didlpi{v,e,t) 

+ / ^m{e - e')pi{v,e,t)pi{v,e' ,t)de' 

^ Jo 

^^de f \m[a{e-e')]pi{v,e,t)p2{v,e',t)de', (?) 

^ Jo 

dp2{v.0.t) ^ 
dt 

— 2u2 [ COS dx -\- sin 6 dy] p2 (r, ^, t) 
^D2dlp2{v,e,t) 

+ / ^m{e - e')p2{v,e,t)p2{v,e' ,t)de' 

^ Jo 

-v^de [ " sm[a{e-e')]p2(v,e,t)pi{v,e',t)de'. (8) 

^ Jo 

Here, a = -2 for Eq. ([5| and a = 1 in case of Eq. Q. 

It can be seen from Eqs. ([l])-(|3|, ([7|), and ^ that 
Uj {j = 1,2) have dimensions of velocity and Dj {j = 
1,2) have dimensions of inverse time. From dimensional 
analysis it follows that the parameters ^i, ^2, and g can 
be measured in units of uf/Di or U2/D2. 

Later in this first study, we will restrict ourselves to 
the case where the particles of the two species feature 
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identical physical behavior. More precisely, this means 
ui =1^2, Di =1^2, and gi = g2' In this case, we can set 
Ui = U2 = 1, 1^1 = D2 = 1, and ^1 = ^2 = 1 without 
loss of generality. The latter can be seen, for example, 
from rescaling time, space, and densities via t t' /Di^ 
r y'ui/Di, Pi p'^Di/gi, and p2 P2^i/^2. As 
a consequence, there are only four different parameters 
remaining, namely 2x2 /'^i, ^2/^1, g/gii and g/g2' For 
the special case of identical physical behavior of the two 
particles there remains only one independent parameter 
that we can define as g' = g/ gi = g/ g2' 



III. LINEAR STABILITY ANALYSIS 

We see from Eqs. ^ and ([8| that pi{r^O,t) = pio and 
P2(r, = p2o is always a solution. Here, pio and p2o 
give the spatially and angularly averaged densities for the 
two species. They are conserved quantities. 

As a first step, we check the linear stability of this triv- 
ial solution with respect to orientational order. This is 
interesting because the points of instability of the trivial 
solution correspond to the onset of collective motion. 

We will not explicitly include the spatial component 
in these considerations. The gradient term only leads to 
an additional imaginary contribution to the listed eigen- 
values. Therefore it does not modify the location of the 
obtained threshold points at which collective motion sets 
in. 



A. Polar and antiparallel alignment 

To check linear stability for the case of polar and 
antiparallel alignment rules between different particle 
species, we must choose a = 1 corresponding to Eq. Q. 
We insert the ansatz 

Pj{r,e,t) = pjo + Pjoe''"'+''\ j = l,2 (9) 

into Eqs. ^ and ([8|. Linearizing in the amplitudes pjo, 
j = 1,2, leads to 

Xpio = - Din'^pio + (^ipio + gp2o)pioSni, (10) 

Ap20 = - D2n^p2Q + (^2P20 + gplQ)p2Q5nl' (H) 

Therefore, only the mode n = 1 can become linearly 
unstable. 

From these equations, we obtain the eigenvalues 

A = ^|(^lplO - Di) + (^2P20 - D2) 
± \/[(^lPlO - ^1) - (^2P20 - D2)f +4^VlOP2o|. 

(12) 

On increasing the average densities, the -eigenvalue 
always becomes unstable first and defines the onset of 
the linear instability. Interestingly, the sign of g does 



not play a role. Therefore, the two vectors of collective 
motion of the two different species can align in a po- 
lar or antiparallel way at onset. Without any coupling 
{9 = 0)7 we correctly recover the behavior of the sepa- 
rate single-component systems 1 and 2 with the respec- 
tive eigenvalues \j = gjPjo — j = 1,2. This case 
leads to the well-known threshold densities p^^ for the 
single-component systems: we find collective motions for 
densities higher than p*Q^ = ^j/dj^ j = 2. 

Clearly, if both pio and p2o are higher than the thresh- 
old single-component densities, we also find collective 
motion in the two-component system. Can collective mo- 
tion also set in, however, if both pio and p20 are smaller 
than the threshold single-component densities? 

Analysis of Eq. (12) shows that for this purpose 



Pio [{gig2 - g^)p2o - ^2^1] < A (^2^20 - ^2) (13) 

is a necessary condition. On the one hand, if g^ > ^1^2, 
P20 can have any value, and collective motion sets in at a 
lower density pio than in the single-component case. The 
threshold density pIq is then given by condition (13). 
On the other hand, if g^ < ^1^2, it follows from Eq. (12) 



that both pio and p2o can be smaller than their threshold 
single-component values and induce collective motion. In 
detail, this can be seen by rewriting Eq. (13) for the case 
that P20 < D2gi/{gig2 - g'^) to 



Pio > 



g2p2Q> 



D2 



gi {g2 - g'^/gi)p\ 



20 



Do 



(14) 



For P20 = we recover the single-component condition 
for pio- When we now increase P207 the (negative) nu- 
merator in Eq. ( 14 ) grows faster than the (negative) de- 
nominator because g2 > g2 — g^ / gi > 0. 

As a result, the two species support each other in start- 
ing to move collectively. In this way, the critical threshold 
densities can be smaller than for the single-component 
systems. Of course, the above analysis analogously ap- 
plies when the two species are switched by changing the 
subscripts 1 2. 



B. Perpendicular alignment 

Considering now a = —2 and again ansatz ([9| in 
Eqs. (0 and ([8|, we obtain 

ApiO = - DiTl^piQ + ^iPloPlO^^nl 

- 2gpiop2oSn2, (15) 

Ap20 = - D2n^p2Q> + ^2P20P20^nl 

- 2gp20pl0^n2' (16) 

Consequently, now the modes n = 1 and n = 2 can be- 
come unstable. 

On the one hand, mode n = 1 becomes unstable if at 
least one of the two averaged densities is higher than the 
threshold one 



o*'^ - ^ 7-12 



(17) 
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These are the same values as for each one-species system 
separately. 

On the other hand, mode n = 2 becomes unstable 
for the density product P10P20 larger than the threshold 
product 



*,2 *,2 

Pio P20 



(18) 



To become unstable before the mode n = 1 and thus 
below the density threshold of the single-component sys- 
tem, we must have p^^ < p*Q^, j = 1,2. For this to be 
possible, the parameters must satisfy the condition 



(19) 



It follows from Eqs. (15) and (16) that the instability 
of the mode n = 2 leads to pio and p20 having opposite 
signs. Consequently, particles of species 1 then on aver- 
age propel perpendicularly to particles of species 2. We 
can say that different species try to evade each other by 
perpendicular alignment of their velocity vectors. It is in- 
teresting to note that the addition of a strongly avoiding 
species can reduce the critical density at which collective 
motion sets in. 



IV. PARTICLE SIMULATIONS 

As we have seen at the end of the last section, the case 
of perpendicular alignment of the velocity vectors of the 
different species is more complex than the case of polar 
and antiparallel alignment. Two alignment modes com- 
pete at onset, corresponding to polar (mode n = 1) and 
nematic (mode n = 2) orientational order. We can iden- 
tify the coupling parameter g as the crucial parameter 
that controls the nature of the orientational mode. To 
get more insight into the corresponding relations, we have 
performed particle simulations of binary self-propelled 
particle mixtures. Here, we mainly report our results 
as a function of the coupling parameter g. According to 
the rescaling procedure pointed out at the end of section 
[m the inter-species coupling parameter g will be given in 
units of the intra-species coupling parameters gi = g2- 

We chose a two-dimensional quadratic simulation area 
of size Lx X Ly with periodic boundary conditions and Nj 
particles of species j {j = 1,2). The quadratic shape was 
used to offer symmetric directions in the case of global 
perpendicular alignment. For simplicity, the velocities 
were set to -ui = -^2 = 1 and the interaction parameters 
between particles of the same species to ^1 = ^2 = 1- 
During each iterating time step, the orientation angle Oi 
of each particle was updated according to Eq. ([2|. In 
contrast to Eq. (|3|, however, particle interactions were 
not completely localized for practical reasons. We chose 
a disk-like environment of radius do around each parti- 
cle, where we mostly used do = 0.05. Interactions were 
considered between particles within this distance. For 
this reason, results from the particle simulations can- 
not be quantitatively compared with the results from the 



Fokker-Planck approach. Apart from that, at each time 
step, we disturbed the orientation angle of each particle 
by an additive Gaussian stochastic noise term according 
to Eq. ([2|. The variance of the Gaussian distribution 
was chosen such that the angular diffusion parameters 
were Di = D2 = I. After that, the advection step was 
performed according to Eq. ([T]) . The time increment was 
set to dt = 0.1. Following our rescaling as mentioned at 
the end of section |TI[ distances such as do are measured 
in units of ui/Di = U2/D2 and time steps in units of 
1/Di = I/D2. 

As an initial condition, a random spatial and angular 
distribution of the particles was chosen. However, spatial 
overlap of the surrounding disks of radius do was avoided 
in the initial configuration. To speed up the calculations, 
a cellindexing method was used [43 . Here, the number 
of cells was adjusted to the size of the simulation area 
and the particle densities. 

To describe the onset of collective motion and the de- 
gree of ordering of the particle velocity vectors, we define 
the usual global order parameters. On the one hand, po- 
lar order is characterized by a vector P j with components 



^^^^ = J^Y.^^'^^^o^ 



Pj.y = ]^IZ^^^^^.' 



(20) 
(21) 



where the sum is over all particles kj of species j (j = 
1,2) and A^2 = N — Ni equals the number of particles 
of species 2. The degree of global polar order is then 
obtained as the magnitude of this vector. 



(22) 



On the other hand, nematic order is described by a trace- 
less symmetric tensor Qj with components 



Qj,xx 
]j,xy — Qj,yx 



-y 



cos 6//e^. 



1 



-y 



(23) 
(24) 
(25) 



which implies Qj^xx = ~Qj,yy (j = 1^2). We obtain the 
degree of global nematic order as 



^j,yy 



(26) 



In general, spatial heterogeneities emerge close to the 
onset of collective motion as noted in the Introduction. 
To first focus on the orientational order itself and re- 
duce the effect of possible spatial heterogeneities, we 
studied relatively small systems. These were of size 
Lx = Ly = 0.83 with Ni = N2 = 117 particles. We 
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FIG. 2. (Color online) Polar and nematic degrees of orient a- 
tional order Pj and Sj as a function of the coupling parameter 
g for preferred polar (g > 0) and antiparallel (g < 0) align- 
ment between the two species j = 1,2. Approximately, the 
two cases are symmetric to each other, with slightly smaller 
values on the antiparallel {g < 0) side for higher values of l^^l- 
The system size was relatively small with Lx — Ly — 0.83 
and Ni — N2 — 117 to reduce the influence of spatial het- 
erogeneities. Results are averaged over 100 independent runs. 
Other parameter values were ui — U2 — 1^ gi — g2 — and 
Di = = 1 during the corresponding particle simulations. 
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FIG. 3. (Color online) Snapshots of particle simulations for 
densities above the onset of collective motion for preferred 
polar (g > 0) and antiparallel {g < 0) alignment between 
different particle species. The system size was Lx = Ly = 10 
with Ni — N2 — 17000 particles, the snapshots were taken 
after 30000 iterations of time step dt = 0.1. Other parameter 
values were ui — U2 — 1^ gi — g2 — and Di — D2 — 1. 
The two species are indicated by red (darker in grayscale) 
and turquoise (brighter in grayscale), and only every 20t/i 
particle position is marked. Clearly, an increasing magnitude 
of coupling 1^1 between the two species decreases the spatial 
heterogeneity in the system. 



varied the coupling parameter g and determined the de- 
grees of orientational order. These were averaged over 
100 independent runs that started from different initial 
conditions. 

To study the occurrence of spatial heterogeneities 
above onset, we turned to larger system sizes. We report 
results for quadratic simulation areas of L^ = Ly = 10 
and Ni=N2 = 17000. 



A. Polar and antiparallel alignment 

The polar and nematic degrees of orientational order 
obtained for the small system sizes are depicted in Fig. |2] 
for each species. Due to the inherent symmetry the two 
species show approximately the same behavior. For each 
particle species separately, the density is above the crit- 
ical threshold density for the onset of collective motion. 
Therefore, we obtain nonzero order parameters already 
for vanishing coupling between the two species at ^ = 0. 

We find polar orientational order within each species 
Pj 7^ that also leads to a nonzero value S'j 7^ 0, j = 1, 2. 
With increasing magnitude of coupling \g\ between the 
two species the orientational order increases. The two 
species support each other in orientational ordering. Ap- 
proximately, the curves for ^ < and ^ > are symmet- 
ric to each other with respect to the line at ^ = 0. At 
higher values of |^| the values for ^ < are slightly lower 
than the ones for ^ > 0. This seems natural since for pre- 
ferred polar alignment the interacting pairs of particles of 
different species move into the same direction. They have 



a longer time of interaction compared to the antiparallel 
case, where they only meet at an instant. 

For larger system sizes, spatial heterogeneities develop 
above the onset of collective motion. Example snapshots 
are shown in Fig. [3] 

We start our discussion with the snapshot for ^ = 0.0. 
In this case the two particle species are decoupled and 
form independent subsystems. Both particle densities 
themselves are above (but not too far above) the criti- 
cal single-particle system density. Therefore, we observe 
the emergence of the traveling stripes or fronts as it was 
found for the single-particle case [30l [31] |33l [36]. The 
stripes move perpendicularly to the direction of their 
elongation with nearly the single-particle speed. It is 
by accident that the stripes for different particles are ori- 
ented almost perpendicularly to each other in this pic- 
ture. 

When we decrease the coupling parameter g to nega- 
tive values, antiparallel alignment between the two par- 
ticle species is preferred. In the snapshot for g = —0.2 
we can see that the stripes are now oriented parallel to 
each other. The two stripes for different species feature 
opposite directions of collective motion. Therefore, from 
time to time they penetrate through each other satisfying 
the antiparallel alignment rule. 

On further decreasing ^, the spatially heterogeneous 
stripe order is destroyed as can be seen from the picture 
for g = —0.4, and the system becomes more spatially ho- 
mogeneous. This process takes place in the regime where 
the degrees of orientational order still strongly increase 
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with 1^1 . 

When we increase g to positive values, we induce po- 
lar alignment between the two particle species. Now the 
stripes for the two particle species tend to move into the 
same direction. As a result, particles of different species 
mix and form one large stripe, as shown by the snap- 
shot for g = 0.6. Again, this compound stripe moves 
perpendicularly to the direction of its elongation. 

Further increasing g destroys the spatial heterogeneity 
also in the case of polar alignment. In the picture of 
^ = 1.0, stripe-like residues are still visible, but they do 
not form a compound object as for the case of g = 0.6. 
For ^ = 5.0 the thin stripes are not visible any more. 

In summary, we found that an increasing magnitude of 
coupling interaction \g\ reduces the spatial heterogeneity. 
We can understand this by analogy to the single-particle 
case. For the latter, it was reported that stripes only 
persist above but close to the onset of collective motion 
[30l [3ll [33] . Far above onset, the single-particle systems 
were found to become spatially homogeneous again. 

Effectively, we observe the same phenomenon in our 
systems when we increase the coupling between the two 
species \g\. It is easiest to see this in the polar case. For 
the value g = 1.0, we effectively obtain a single-species 
system of twice the density as for g = 0.0. Therefore, we 
have effectively increased the particle density to far above 
its threshold value for the onset of collective motion. This 
significantly reduces the spatially heterogeneous ordering 
into stripes. The latter becomes obvious when we com- 
pare the snapshot for ^ = 1.0 to the cases of ^ = 0.0 and 
g = 0.6. 

Apparently, the effect is much stronger in the case of 
preferred antiparallel alignment between different parti- 
cle species when compared to polar alignment. This is ex- 
pected for the following reason. In the polar case, where 
the two species form one compound stripe, the interac- 
tion between particles can be effectively reduced. The 
latter is possible by an increased elongation of the stripe 
or by splitting into several stripes. Since all of these 
objects move into the same direction, they do not meet 
each other and consequently the frequency of interaction 
is lower. In the antiparallel case, the stripes for different 
species move into antiparallel directions and from time 
to time must penetrate through each other. When these 
events occur, the density of particles interacting is dou- 
bled within the stripes and therefore remains high above 
the threshold density for the onset of collective motion. 



B. Perpendicular alignment 

Again, we first focus on the small systems, where the 
role of spatial heterogeneities is reduced. As expected 
from section [lll[ this case of perpendicular alignment is 
qualitatively different and richer in phenomena than the 
corresponding parallel case. The results for the degrees 
of polar and nematic order are depicted in Fig. |4] We 
can see that the two particle species show quantitatively 




FIG. 4. (Color online) Polar and nematic degrees of orient a- 
tional order Pj and Sj as a function of the coupling parame- 
ter g for preferred perpendicular alignment between the two 
species j = 1, 2. The values Sj^caic were calculated solely from 
the magnitude of the corresponding Pj values. A transition 
from predominantly polar to non-polar nematic order is ob- 
vious around g = 2. Technical details as given by the caption 
of Fig. [I 



the same behavior, as expected for reasons of symmetry. 

For small values of g polar order dominates within each 
species. Here, the magnitudes of the degrees of polar 
orientational order Pj {j = 1,2) are close to 1. As we 
can see from Fig. [4j the values of Pj slightly increase 
with increasing g. Apparently the scattering between 
particles of different species enhances the polar ordering 
within each species for small values of ^. In addition, 
we could show that the nonzero degree of nematic order 
in this regime is only due to polar orientational order. 
For this purpose, we assumed a Gaussian distribution 
of the orientational angles Ok^ around their mean value 
for each species. This assumption is corroborated by the 
observations from the simulations. For each value of g, 
the Gaussian distribution is completely determined via 
the corresponding value of Pj . We then use this Gaussian 
angular distribution to calculate the degree of nematic 
orientational order, which we call Sj^caic As depicted in 
Fig. [4] the values thus obtained are identical with those 
of Sj that were extracted directly from the simulations. 

Around the value g = 2 polar orientational order 
breaks down. Eq. (19) predicts a threshold value of g 
at which the mode n = 2 can become unstable first. 
Indeed, we observe that the degrees of nematic orien- 
tational order Sj further increase with increasing values 
of ^ > 2. The small dip in Sj reflects the drop in polar 
orientational order, which is then compensated by truly 
nematic (mode n = 2) orientational order. 

For the larger system sizes that feature spatial hetero- 
geneities in the form of stripe textures we depict example 
snapshots in Fig. [Sj The case is similar to the one of pre- 
ferred antiparallel alignment that was shown in Fig. [3] for 
g>0. 

The state of decoupled species at ^ = 0.0 is identical to 
the one in Fig.jS] Already for g = 0.2 the spatial ordering 
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FIG. 5. (Color online) Snapshots of particle simulations for 
densities above the onset of collective motion for preferred 
perpendicular alignment between different particle species. 
Technical details are the same as given by the caption of 
Fig. [3] Again, an increasing magnitude of coupling g between 
the two species decreases spatial heterogeneity in the system. 



into stripes is noticeably reduced. The stripe objects 
move approximately perpendicularly to each other. For 
g = 0.4 the stripes have dissolved. 

As in the antiparallel case, the three snapshots fall into 
the regime where the overall degrees of orientational or- 
der still strongly increase with the magnitude of g. In 
contrast to the antiparallel case, the stripes never com- 
pletely overlap along their elongation. However, they are 
always in contact at one crossing intersection. There, 
the density of interacting particles is again doubled and 
therefore far above the onset value for collective motion. 



V. NUMERICAL SOLUTIONS OF THE 
FOKKER-PLANCK EQUATIONS 

The results of the particle simulations compare well to 
the numerical results obtained from the continuum equa- 
tions. To solve Eqs. ^ and Q numerically, we used 
a finite differencing scheme. In analogy to the particle 
simulations, the spatial calculation grid was quadratic of 
size Nx X Ny with periodic boundary conditions. is 
the number of possible angular orientations considered. 
It was chosen such that the convolution integrals in the 
angular distributions could be efficiently evaluated via 
fast Fourier transforms [44]. We found that a first order 
upwind scheme for the first order spatial derivatives re- 
produces well the physical properties of the system. The 
time step must be small enough to conserve the overall 
particle densities. 

As for the particle simulations, we focus on two cases. 
First, we investigate spatially homogeneous solutions by 
setting Nx = Ny = 1 {N^ = 128). After that, the influ- 
ence of spatial degrees of freedom is taken into account, 
where we mostly used Nx = Ny = 32 {N^ = 32). The 
parameters were set as for the particle simulations, ex- 
cept for the velocities {ui = U2 = 0.1, ^1 = ^2 = 1, 
Di = D2 = 1), and we varied the interaction param- 
eter g. For our choice of parameter values, we obtain 
from Eq. (17) a critical single-species system density of 



p*Q = 1 {j = 1,2). We report results for three charac- 
teristic scenarios: (a) both values of the mean densities 



Pjo {j = 1,2) are above the critical single-species system 
densities (pio = P20 = 1-5); (b) one of the mean densities 
is above the critical single-species system density and one 
is below (we consider pio = 1.5, p2o = 0.5 for spatially 
homogeneous solutions and pio = 1.1, P20 = 0.5 for spa- 
tially inhomogeneous ones); and (c) both mean densities 
are below the single-species system density (pio = P20 — 
0.5). For each value of ^, we initialized the densities on 
the Nx X Ny X N^ sized calculation grid by the mean 
densities plus a random number of Gaussian distribution 
and small amplitude. 

In the continuum picture, we obtain the local orien- 
tational order parameters by taking the moments of the 
densities pi(r, 6>,t) and p2(v^0^t) with respect to the an- 
gular distributions, 

Cj{v,t) = / pj{r,e,t)dd, (27) 
Jo 

cj{r,t)P,{v,t) = j^J ""^^l ^ p,ir,0,t)d0, (28) 

c,(r,t)Q,(r,t) = y^ (.cos^sin^ sin'O - I 

X pj{r,0,t)dO, (29) 

j = 1,2. Here, Cj(r, t) gives the local particle number 
density. Pj(r,t) corresponds to the /oca/ polar alignment 
vector, whereas Qj(r, t) is the /oca/ nematic order param- 
eter tensor for each species j = 1,2. At each position r, 
the local degrees of orientational order follow in analogy 
to Eqs. (22) and (26). To obtain global degrees of orienta- 
tional order, we took the spatial averages over the system 
size. Discretized versions of these definitions were used 
for the numerical implementation. 



A. Polar and antiparallel alignment 

First, we confine ourselves to spatially homogeneous 
solutions. We show characteristic results for preferred an- 
tiparallel alignment between particles of different species 
in Fig. [g] as a function of |^| (results for polar alignment 
follow approximately analogously for ^-values of oppo- 
site sign, compare also Fig. [2|. The degree of nematic 
order is nonzero due to polar orientational order and not 
depicted in the figure. 

For pio = P20 = 1.5 both species densities are above 
the critical one-component density. There is collective 
motion already without coupling at ^ = 0. The relative 
angular orientation between the two species is not yet 
fixed, however. For nonzero values of |^|, the polar or 
antiparallel orientation of the two species velocities sets 
in. The two species support each other in orientational 
ordering with increasing values of g. 

The asymmetric case, where pio = 1.5 is above and 
P20 = 0.5 is below the threshold density, is interesting 
for the following reason. For both species, orientational 
order increases due to the coupling to the other species 
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FIG. 6. (Color online) Polar degrees of orientational order 
Pj for the two species j = 1, 2 as a function of the coupling 
parameter g. The results were obtained from spatially homo- 
geneous numerical solutions of the Fokker-Planck equations 
for preferred antiparallel orientational order between the two 
species. Different mean particle densities were considered. 
Other parameter values were ui = U2 = 0.1, gi = g2 — 1, and 
Di=D2 = 1. 



with increasing |^|-values. For the given density values, 
P2 starts to grow already at ^ = 0. Remarkably, the de- 
gree of orientational order of the more diluted species P2 
even exceeds the one of the denser species Pi for strong 
coupling between the two species. Surprisingly, for high 
values of |^|, it is P2 that asymptotically approaches the 
degree of order that was reached for pio = p2o = 1.5. In 
contrary. Pi asymptotically approaches the lower degree 
of order reached for pio = p20 = 0.5 (see below), de- 
spite the higher mean density pio = 3p20- So Pi and P2 
behave oppositely to what would be expected from the 
corresponding mean densities. However, we note that, 
at high values of |^|, orientational order of one species 
is predominantly achieved by interactions with the other 
species. In this way, it is the density of the other species 
that determines the asymptotic degree of ordering. 

In the case of pio = p2o = 0.5 there is no collective 
motion for ^ = 0. Only at \g\ = 1 collective motion sets 
in for both species simultaneously. At this value, the two 
species form an effectively single-component system, so 
that the densities add up to an effective density pio + 
P20 = 1. This value is the threshold density for the onset 
of collective motion in a corresponding single-component 
system [compare, e.g., Eq. (17)]. 

Second, we investigated the spatial heterogeneities that 
appear for the larger system sizes. Qualitatively, our re- 
sults obtained from the particle simulations in the pre- 
vious section are confirmed. They correspond to the 
case pio = p2o = 1-5, in which both mean densities are 
above the single-species threshold density. Spatial het- 
erogeneities in the form of stripes appear in the density 
profiles. An example is shown in Fig. [7](a) and (b). 

Again, we find that an increasing magnitude of the 
coupling parameter \g\ dissolves the stripes. As for the 



FIG. 7. (Color online) Spatial density distributions obtained 
from numerical solutions of the Fokker-Planck equations for 
mean densities pio = p2o = 1.5. Upper panels [(a) and (c)] 
correspond to species 1, lower ones [(b) and (d)] to species 2. 
Left panels [(a) and (b)] follow from preferred polar alignment 
between the different species with g = OA and show parallel 
orientation of the resulting stripes. Right panels [(c) and 
(d)] follow from preferred perpendicular alignment with g = 
0.6 and feature perpendicular orientation. Other parameter 
values were ui = U2 = 0.1, gi =^2 = 1, and Di = D2 = 1. 
The numerical grid size was Nx x Ny = 32^ lattice points of 
distance dx = 0.5 with N^^ = 32 angular orientations each, 
and the equations were iterated 3 x 10^ times with step size 
dt = 0.001. Brightness increases with density and has been 
rescaled to maximize the spatial density contrast. 



particle simulations, the necessary values for \g\ are much 
smaller for antiparallel than for polar alignment between 
different species. The numerical solutions of the Fokker- 
Planck equations offer a simple method to determine this 
value of 1^1 : the difference between the largest and the 
smallest density value decays to zero when \g\ destroys 
the spatial heterogeneity. We find a value of |^| ~ 4.8 
in the polar and \g\ < 0.2 in the antiparallel case. In 
contrast to the particle simulations, we often observe the 
velocity vectors to align along the stripe direction for an- 
tiparallel interactions. The density profiles of the two 
stripes of different species then are stationary and over- 
lap. This increases the interaction time between the two 
species densities. 

The asymmetric case of pio = 1.1 and p2o =0.5 shows 
interesting phenomena. For the majority species 1 the 
mean density is above onset. Consequently collective mo- 
tion sets in and spatial heterogeneities occur. Through 
the interaction of strength ^, collective motion can also 
be induced in the minority species 2. However, this hap- 
pens only at positions where the density Ci(r, t) is high 
enough. 

For polar alignment interactions. Fig. [s] (e) shows the 
final density profile of a stripe of species 1. Within the 
stripe region, Ci(r, t) is high and induces collective mo- 
tion in species 2. At these locations, the material of both 
species propels into the same direction. As a result, spots 
of high density C2(r,t) follow the ones of high density 
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FIG. 8. (Color online) Evolution of the spatial density dis- 
tribution obtained from numerical solutions of the Fokker- 
Planck equations for mean densities pio = 1.1 and p2o = 0.5 
in the case of polar alignment. The coupling strength be- 
tween the two species was g = 2.2. Upper panels [(a), (c), 
and (e)] correspond to species 1, lower ones [(b), (d), and (f)] 
to species 2. Shown snapshots were taken at the following 
numerical times: (a) and (b) 10; (c) and (d) 50; (e) and (f) 
3000. At each depicted instant the density map for species 2 
is a copy of the map for species 1. Further technical details 
as given by the caption of Fig. [7| 



ci(r, t). In this way, the density map of C2(r, t) depicted 
in Fig. [s] (f) becomes a copy of the one of ci(r, t). Con- 
sequently, an induced overlapping stripe of species 2 is 
generated with polar alignment of the collective velocity 
vector. These statements even hold at early times of the 
ordering process when the stripe textures have not yet 
developed. An example is given by the time series in 
Fig.|(aHf). 

It is interesting to note, however, that for the antipar- 
allel case the density profiles are inverted. This is illus- 
trated by Fig. [9] (a) and (b), where a stripe of species 1 
is moving to the left. The motion to the left leads to the 
sharp front and the fuzzy tail of the stripe in panel (a). 
In the region of the stripe, the density Ci(r, t) is high so 
that it can induce collective motion in species 2. This 
happens via the antiparallel alignment interaction ^ < 0. 
Consequently, the vector of collective motion of species 
2 points into the opposite direction, i.e. to the right. In 
this way, material of species 2 at the head of the mov- 
ing stripe is "pumped" through the stripe of species 1. 
Behind and outside the stripe of species 1, the density 
ci(r, t) is so low that it cannot induce effective collec- 
tive motion in species 2 any more. Thus the material of 
species 2 is not advected once it has passed the stripe. 
It gathers behind the stripe area. In effect, this leads to 
the inverted stripe density profile shown in Fig. |9](b). 

As in the spatially homogeneous case, collective motion 
in the system pio = P20 = 0.5 sets in at values \g\ > 1. 
Stripes develop above onset. 



FIG. 9. (Color online) Spatial density distributions obtained 
from numerical solutions of the Fokker-Planck equations for 
mean densities pio = 1.1 and p2o = 0.5. Upper panels [(a) and 
(c)] correspond to species 1, lower ones [(b) and (d)] to species 
2. Left panels [(a) and (b)] follow from preferred antiparallel 
alignment between the different species with g — —0.2. The 
stripe in panel (a) travels to the left and features the typical 
sharp front and fuzzy tail. Material of species 2 is advected 
through the stripe in opposite direction. Right panels [(c) 
and (d)] follow from preferred perpendicular alignment with 
g = 2.2. Again the flow field in panel (c) is oriented to the left. 
Material of species 2 is expelled from the stripe region to the 
top and bottom. Both cases lead to inverted spatial density 
profiles for species 2. Further technical details as given by the 
caption of Fig. [7| 



B. Perpendicular alignment 

For preferred perpendicular alignment between parti- 
cles of different species we again focus on spatially ho- 
mogeneous solutions first. Typical results are depicted 
in Fig. 10 We only show the degree of nematic order Sj 
Degrees of polar orientational order remain 



12] 



(i 

negligibly small, Pj ^ 0, for corresponding mean den- 
sities pjo = 0.5. If the mean density is poj = 1.5, the 
Pj curve has qualitatively the same shape as the ones in 
Fig. |4j with the strong descent located at the dip of the 



respective curve for Sj in Fig. 10 

At mean density values pio = P20 = 1.5 collective mo- 
tion of polar order dominates at low coupling strength g. 
This polar order breaks down where we find the dips in 
the curves of Sj (j = 1, 2) in Fig. flOj In contrast to the 
particle picture, where Pj ^0.1 after this transition (see 
Fig.[4|, we now find values close to zero. At higher values 
of ^, truly nematic order dominates. Our linear stability 



analysis from section [III] does not provide a quantitative 
measure for the location of the transition because of the 
nonzero amplitudes of the polar order parameters below 
the transition. Differences when compared to the par- 
ticle simulations illustrate the idealized character of the 
mean field continuum approach, where finite interaction 
radii and free paths between particle interactions were 
not taken into account. 
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FIG. 10. (Color online) Nematic degrees of orientational or- 
der Sj for the two species j = 1, 2 as a function of the coupling 
parameter g. The results were obtained from spatially homo- 
geneous numerical solutions of the Fokker-Planck equations 
for preferred perpendicular orientational order between the 
two species. Different mean particle densities were considered. 
Other parameter values were ui = U2 = 0.1, gi = g2 — 1, and 
Di=D2 = 1. 



The asymmetric case of pio = 1.5 and p2o = 0.5 fea- 
tures the same interesting effect as in the previous sub- 
section for antiparallel alignment. Species 1 shows col- 
lective motion already at zero coupling ^ = 0. Polar ori- 
entational order dominates and increases with increasing 
coupling strength. The polar orientational order breaks 
down at the location of the dip in the Si curve in Fig.[lO 
After that, truly nematic order prevails. Si then asymp- 
totically approaches the curves of Sj for pio = p2o = 0.5 
(see below) although here pio = 1.5. In contrast, orien- 
tational order for species 2 increases from zero with in- 
creasing coupling strength. Polar order never plays a sig- 
nificant role. The value of 5*2 asymptotically approaches 
the value of Sj for pio = P20 = 1.5, despite the fact that 



P20 = 0.5. In analogy to the previous subsection, here 
Si and 5*2 behave oppositely to what would be expected 
from the corresponding mean densities. Again, it is the 
density of the other species that determines the asymp- 
totic degree of ordering of one of the two species. This 
is because at high values of g orientational order of one 
species is predominantly achieved by interactions with 
the other species, and not by interactions with particles 
of the same species. 

In the case of pio = P20 = 0.5, collective motion sets 
in for both species simultaneously at ^ = 4. This value 
is predicted by the linear stability analysis via Eq. (18). 



Above this threshold, truly nematic orientational order 
dominates. 

When we turn to the larger system sizes, spatial inho- 
mogeneities can arise. For pio = p2o = 1.5, both mean 
species densities are above the critical single-species den- 
sity. Qualitatively we then find similar behavior as for 
the particle simulations in Fig. |5j Stripes develop that 
are oriented perpendicular to each other as illustrated by 



the example in Fig.[7|(c) and (d). Similarly to the case of 
antiparallel alignment, the velocity vectors of collective 
motion were mainly oriented along the stripe direction. 
Again, increasing the coupling strength between the two 
species dissolves the stripes. We found that the system 
turns spatially homogeneous for ^ > 1. 

For pio = 1.1 and p2o = 0.5 we observe the same ef- 
fect as for the antiparallel alignment in Fig. [9] Since the 
density for species 1 is above the critical density, collec- 
tive motion sets in and spatial heterogeneity in the form 
of a stripe develops. We find these stripes for g < 4.8. 
At spots of high density ci(r, t), i.e. within the stripes, 
collective motion is induced in species 2 due to the inter- 
action of strength g. An example result is illustrated in 
Fig.[9](c) and (d). In this situation, the velocity vector of 
species 1 is oriented along the stripe. Therefore, material 
of species 2 is pumped through and out of the stripe. The 
consequence is a depletion of material of species 2 within 
the stripe region and again an inverted density profile for 
species 2, as shown by Fig. |9](d). 

Finally, when pio = P20 = 0.5, we recover the critical 
value of ^ = 4 that was already found for the spatially 
homogeneous case in Fig. 10 Above this value, and 



5*2 become nonzero and collective motion develops. Inter- 
estingly, we found that the systems do not turn spatially 
heterogeneous directly above onset. At the same time, no 
polar orientational order was detected. Only for values 
g > 4.8 we observed stable spatial heterogeneities. The 
emergence of these spatial heterogeneities, however, was 
coupled to the development of polar orientational order 
within each species. This observation is in contrast to 
the spatially homogeneous solution, where in this regime 
we found Pi ^ ^ P2. It seems that for the case of 
perpendicular alignment and pio = p2o = 0.5 spatial 
heterogeneities are coupled to the evolution of nonzero 
values of Pi and P2. 



VI. MACROSCOPIC CONTINUUM 
EQUATIONS 

In this section, we derive macroscopic hydrodynamic- 
like equations from the Fokker-Planck equations ^ and 
([8| and briefly discuss their regime of validity. We ob- 
tain the characteristic macroscopic variables by taking 
the moments of the densities pi(r,^,t) and p2(r, ^,t) 
with respect to the angular distributions as given by 
Eqs. (27)-(29). Assuming one single mass rnj for par- 
ticles 01 eacETspecies, Cj(r,t) is proportional to the mass 
density. Pj(r,t) gives the local polar alignment vector, 
whereas Qj (r, t) characterizes the local nematic order for 
each species j = 1,2. Since the magnitude of the ve- 
locity is fixed for each species, Cj(r, t)Pj(r, t) is propor- 
tional to the momentum density mjCj(r, t)vj(r, t). Here, 
the macroscopic velocity field Vj(r,t) can be obtained by 
averaging over all velocity vectors of particles located at 
time tin a surface element at position r. The correspond- 
ing particle velocity vectors are given by Eq. g. This 



12 



leads us to the expressions 



Cj(r, t)vj(r, t) = 2u^ 



Jo 



l-fo ] P,ir,e,t)d0, (30) 



where j = 1,2. 

The dynamic equations are derived by taking the mo- 
ments of Eqs.J?) and ([8| as given by the right-hand sides 
of Eqs. (27)-(29). Higher order angular moments are not 
taken into account. We thus expect quantitative devi- 
ations of the results obtained from the hydrodynamic- 
like equations when compared to direct solutions of the 
Fokker-Planck equations. This becomes more and more 
severe when densities are significantly higher than the 
threshold values. The equations that we list contain 
terms up to second order in the particle densities ci and 

In both cases of alignment (a = 1 and a = —2), the 
zeroth moment of Eqs. ^ and (|8| leads to the continuity 
equation for each species, 



dCj{T,t) 



,(r,t)v,(r,t)] 



dt ~ ' L-n^^y^.v 

= -2^i,V-[c,(r,t)P,(r,t)], (31) 

j = 1,2. We find qualitative differences for the higher 
order angular moment equations for the different cases 
of alignment. 



A. Polar and antiparallel alignment 

In the case of a= 1, Eq. ^, we derived the following; 
macroscopic equations. We always refer to local particle 
densities as well as local polar and nematic alignment 
order parameters. For brevity, however, temporal and 
spatial coordinates (r, t) are not explicitly noted. 

The temporal evolution of the polar alignment vector 
Pi is given by 



S,(ciPi) 



2ui 
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-P2-P2 Qi 



from which the equation for P2 follows from switching 
all subscripts 



St(c2P2) ~ 1 ^ 2 



(33) 



in Eq. (32). 



Similarly, the temporal evolution of the nematic order 



parameter Qi is obtained as 

Ui 



dt{ciQi) 



[V-(ciPi)]I 



V(ciPi) + {V(ciPi)}' 



2 

4I)iCiQi 

^cl [2PiPi-P?l] 



27r 



C1C2 [PiP2+P2Pl-(Pl-P2)I]. 



(34) 



Here, ^ denotes the transpose of the superscripted gradi- 
ent matrix, and I corresponds to the unity matrix. The 
equation for Q2 follows again from switching all sub- 
scripts 



9t(c2Q2) ~ 1 ^ 2 



(35) 



in Eq. (|34). 

We can see that diffusion tends to reduce both po- 
lar and nematic order through the terms proportional to 
Dj > 0, j = 1,2. This term is obtained from the second 
partial derivative dg in Eqs. ^ and Consequently, 
it grows quadratically in the angular order of the consid- 
ered mode. This allows to approximately neglect higher 
modes of orientational order as mentioned above. 

These equations are now analyzed with respect to the 
threshold for the onset of collective motion. For the rea- 
sons noted in section 



HI 



we restrict ourselves to the 
spatially homogeneous case. I.e. we neglect the gradient 
terms in Eqs. (32)- (35). Close to threshold, the second 
angular momenta, represented by Qi and Q2, relax faster 
than the first angular momenta Pi and P2. Therefore, 
we set the partial time derivatives on the left-hand sides 
of Eqs. (34) and (35) to zero and solve for the stationary 
spatially homogeneous values of Qi and Q2, 



Q^'--4k{i^^ [2PiP.-Pfl] 



+ ^C2 [PiP2+P2Pl-(Pl-P2)I] 



and 



Q 



2, St 



(36) 



(37) 



(32) Inserting into Eqs. (32) and (33) leads to 



at(ciPi) « Ci (I^Ci - Di) Pi + |^CiC2P2 



11 2 

Ci (^iCiPi +^C2P2) Pi 



and 



4Di 27r2 



9t(c2P2) ^1^2. 



(38) 

(39) 
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From a linear stability analysis of these equations around 
the non-ordered state Pi=P2 = 0, we obtain the same 
eigenvalues as those in Eq. (12). This leads to the same 



threshold values for the onset of collective motion as de- 
rived in section ITITl 

Crossing the threshold for the onset of collective mo- 
tion, it is important to note that the system of Eqs. (38) 



and (39) is not stable for large values of the average 



densities. This is already the case for the spatially ho- 
mogeneous single-species scenario that we obtain from 
Eqs. (38) and (39) by setting, for example, C2 = 0. For a 



single-species system, the analog to Eq. (38) reads 



5i(ciPi)«Ci(|^Ci-Z)i 
1 



Pi 



1 



4Di 27r2 



(40) 



Despite the stabilizing cubic term in Pi, a systematic 
linear stability analysis shows that the static solution be- 
comes linearly unstable at density values 



ci > Gtt 



91 ' 



(41) 



Also numerical solutions were found to diverge beyond 
this value. A further analysis demonstrates that den- 
sities satisfying Eq. ( [il] ) lead to the relation Si > Pi. 
This appears unphysical in the case of polar alignment 
between the particles of a single species. 

In the two-species case, as expected. Pi and P2 show 
polar alignment for ^ > and antiparallel alignment for 
^ < 0. The general expressions for the static solutions of 
Eqs. ( 38 ) and ( 39 ) are very lengthy, so we do not list them 



here. Instead, in the following, we confine ourselves to the 
special symmetric case of identical species that interact 
with each other, i.e. Ci = C2, ^1 = ^2, and Di = D2. 

Assuming Pi = P2 and ^i = 5*2 for the degrees of polar 
and nematic orientational order, respectively, we obtain 
the trivial solution Pi = P2 = 0, or 



{gi + \g\)ci - 2TrDi 



The latter implies 

Si = 82 = I 



(gi + \9\rcl 



2ttDi 

{gi + \g\)ci' 



(42) 



(43) 



As we can see from Eq. (42 ), the nontrivial solution exists 
for densities 



ci 



C2 > 



27rP>i 

9i^\9\' 



(44) 



We performed a linear stability analysis in Pi = P2 
and ^1 =5*2, in the angular orientations of Pi and P2, 
as well as in the angular orientations of the principal axes 
of Qi and Q2. This linear stability analysis confirmed the 
critical density values of Eq. ( 44 ) . They are the analog to 



the critical particle densities for the single-species system 



as they follow from Eq. ([T7|). Below this density value, 
the system disorders to the state Pi = P2 =0. On the 
other hand, we found that the solution becomes linearly 
unstable for densities 



ci = C2 > Qn- 



Di 



91 + 1^1 



(45) 



This is the analog to expression (41) of the single-species 



case, where the macroscopic equations are not stable any 
more, gi is replaced by the stronger coupling gi + \g\. 



B. Perpendicular alignment 

The situation can be manifestly different in the case of 
perpendicular alignment, a = — 2, of Eq. ([5|. Following 
the same procedure as in the previous section, we now 
find 



dticiPi) 



2u^ 



V-(ciQi) + ivci 



AciPi 



91 2 

TT 

9 



-Pi -Pi Qi 



C1C2P1 • Q2 



and 



dt{c2P2) ^1^2 

for the polar alignment vectors as well as 

Ui 



(46) 



(47) 



dt{ciQi) 



[V-(ciPi)]I 



V(ciPi) + {V(ciPi)r 



Ui 

2 



gc? [2PiPi-Pi^l] 
4P>iCiQi - -C1C1Q2 

TT 



and 



at(c2Q2) 



(48) 



(49) 



for the nematic order parameters. 

In contrast to the previous case of polar and antiparal- 
lel alignment, the nematic order parameters Qi and Q2 
are now explicitly coupled in Eqs. (|48|) and (49). This is 



a consequence of the different angular interaction poten- 
tial ([5| between the two species. It is of second order in 
the angular momenta. 

Since the angular orientation of one of the polar vector 
order parameters is arbitrary, we can reduce the number 
of independent variables to seven. Still, however, a gen- 
eral systematic stability analysis is out of reach and we 
discuss special stationary solutions in the following. 

Proceeding in the same way as in the previous subsec- 
tion, the stationary spatially homogeneous values of Qi 
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and Q2 now become 



167r^Z)iZ)2 — g'^ciC2 
X I251D2C1 [2P1P1 -P?!] 



992 2 



[2P2P2 



and 



Q 



2, St 



1 ^ 2- 



(50) 



(51) 



We can see that the expression diverges for C1C2 = 
167T^DiD2/ . This indicates that another solution sets 
in at such density values. Indeed, from Eq. ([l8|, our lin- 
ear stability analysis shows us that the second mode of 
orientational order becomes unstable at these densities. 
It corresponds to purely nematic order. For the moment, 
we confine ourselves to the density regime below this di- 
vergence. 



C1C2 < 



167r'^DiD2 



(52) 



First, we assume that one of the two densities is so 
small that only the other species moves collectively. We 
choose Pi 7^ = P2. For the magnitude of polar orien- 
tational order we find 



(51 ci - 27rDi)(167r2DiD2 



5^CiC2) 



gicf{'iTrD2gi -§"^02) 



(53) 



If all terms on the right-hand side are positive, the so- 
lution exists. This is the case for ci above the critical 



single-species density and according to Ineq. (52) if the 
denominator is positive. The latter is true if 



C2 < 



4.7rD2gi 



(54) 



Otherwise Pi diverges and then becomes imaginary. 

The corresponding degree of nematic orientational or- 
der reads 



Si = 



7rZ?2(giCi -27rDi) 
Cl{'iTTgiD2 - g^C2)' 



(55) 



with the principal axis of Qi oriented parallel to Pi. 
Interestingly, the coupling induces nonzero nematic ori- 
entational order of species 2 that is of strength 



^2 = 



_g{giCi^iiDi) 
4{47TgiD2 



■^^C2) 



(56) 



The principal axis of Q2 is oriented perpendicular to the 
one of Qi and to Pi. 

Looking for solutions of non-vanishing polar order for 
both species. Pi 7^ 7^ P2, we find that the relative 



angle between the two vectors is ± | . The magnitudes of 
the vectors are given by 



1 



giciigm - g'^) 

X < (47r^2^i - ^^ci)(^ici - 27rDi) 

-^(4^I)i-^iCi)(^2C2-2^I)2)| (57) 



and 



1 ^ 2. 



(58) 



This solution trivially exists, if the signs of all the terms 
in brackets are positive. Above the single-species thresh- 
olds Cj > 27TDj/gj {j = 1,2), this is achieved if simulta- 
neously ci < 4:7Tg2Di/ g'^ ^ C2 < 4:7TgiD2/ g^ [same condi- 
tion as Ineq. (54)], 



and 



9j 



g < gig2- 



1,2, 



(59) 



(60) 



Indeed we numerically found divergence of the set of 
Eqs. (46)-(49) when the latter condition was violated. 



The corresponding nematic order parameters read 



^1 



^2(^1^1 - 27rDi) + ^(^2^2 - 27ri:>2) 

ci{gig2-g'^) 



and 



1 ^ 2, 



(61) 



(62) 



where we assumed that the principal axis of Qj is parallel 
to Pj to choose the sign of Sj {j = 1,2). 

At the threshold indicated by Eq. (52), namely C1C2 = 
167r^DiD2/g^ ^ the spatially homogeneous part of equa- 
tions ([46|)-(49) has the stationary solution Pi = = P2 
as well as arbitrary Qi and Q2. In order to make a 
statement about the solution above the threshold, non- 
linear terms in Qi and Q2 are necessary. We therefore 
rederived the set of equations (46)-(49) including orien- 



tational moments up to fourth order. The third and the 
fourth modes were used to close the equations. It turned 
out that in the case that is relevant here. Pi = = P2 



namely, Eq. (48) is supplemented by an expression 



47r2Di 



c^^(Q2:Q2)Qi 



(63) 



and Eq. (49) accordingly by a corresponding expression 
of 1 ^ 2- From this it can be shown that a stationary 
spatially homogeneous solution of truly nematic order 
Pi = = P2 and Qi 7^ 7^ Q2 exists above the thresh- 
old given by Eq. (52) as expected. 
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VII. CONCLUSIONS 

In this paper, we have studied the case of binary mix- 
tures of self-propelled particles. We started from a min- 
imal model in which the magnitude of the particle ve- 
locities is kept constant (see, e.g., ref. [25 for the single- 
particle case). The orientational order between the par- 
ticles is adjusted through local pairwise interactions. Po- 
lar orientational order is preferred for particles of the 
same species. For particles of different species, we have 
investigated the cases of preferred polar, nematic, and 
perpendicular alignment interactions. 

We started from the Langevin equations in the parti- 
cle picture. From these, we derived mean field continuum 
equations of the Fokker-Planck type. The onset of collec- 
tive motion and the nature of corresponding solutions in 
the binary self-propelled particle mixtures were studied 
through a linear stability analysis and numerical inves- 
tigations of the particle and continuum equations. Fur- 
thermore, we derived macroscopic continuum equations 
and analyzed corresponding stability ranges. 

It turned out that the interaction between the two 
species can reduce the threshold densities for the onset of 
collective motion to values below the ones for the single- 
species case. If one of the two species has a density above 
this threshold, it moves collectively and can induce collec- 
tive motion also in the other species even if the latter has 
a density below the single-species threshold. In the case 
that both densities are below the single-species thresh- 
old, interaction between particles of different species can 
nevertheless induce collective motion within each species 
for all three alignment rules investigated. 

Above, but close to the onset of collective motion, spa- 
tial heterogeneities in the form of stripe-like fiocks emerge 
in the density profiles. In the most interesting case, one 
of the species densities is above the single-species thresh- 
old and the other is below. Then the first species devel- 
ops collective motion and spatial heterogeneities. Where 
its density is high, it can induce collective motion in the 
second species. Depending on the alignment rules be- 
tween different species, this can lead to identical or in- 
verted density maps for the two species. Increasing the 
coupling strength between the two species dissolves the 
spatial inhomogeneities. 

For the case of preferred perpendicular alignment, we 
find a competition between polar and truly nematic order 
as a function of the strength of orientational coupling 
between the two species. This competition also infiuences 
the development of the spatial inhomogeneities. 

When we are looking for the connection of our study to 
the experimental investigation of real systems, we have 
to keep in mind that our results were obtained for two 
spatial dimensions. We should therefore confine ourselves 
to at least quasi two dimensional systems. Candidates 
for the latter are thin films of motile bacteria colonies 
(monolayers in the ideal case) at air-water surfaces or on 
substrates. For such cultures formed by Bacillus subtilis 
it has been shown that only a fraction of the cells is 



motile, the other cells are non- motile [45] [46]. The size 
of this fraction is controlled by the genetic location of the 
gene responsible for the production of a certain protein 
[45]. In our model system, this situation corresponds to 
the limiting case in which both species are formed by 
the same bacterium, one of the species featuring zero 
motility, the other propelling with non-zero velocity. 

Natural systems in which the latter situation is often 
observed are bacterial biofilms. These are communities 
of microorganisms attached to surfaces. Also biofilms 
of Bacillus subtilis were demonstrated to feature cellu- 
lar differentiation so that only part of the cells are motile 
[45l|47] . During the biofilm development, a fraction of the 
initially motile cells starts to take over a different task 
and forms non-motile sub-communities [4T. Although 
biofilms are typically extended on a surface, their thick- 
ness cannot be neglected and can even feature a spatial 
organization in different layers [47 . However, very thin 
films can be produced for example by letting a biofilm of 
motile bacteria grow in upstream direction in fiow cells 

m- 

Despite their abundance in nature and their much 
higher clinical relevance, multi-species biofilms have only 
recently been moved into the focus of investigation [49] , 
In this case, different species within the biofilm can in- 
teract via quorum sensing and metabolic cooperation or 
competition. Depending on synergistic or antagonistic 
interactions between motile bacterial species, different 
rules of alignment may result when their swarming be- 
havior is investigated. One step into this direction was 
performed by a study on a community of two different 
species of motile bacteria, in which the authors also fo- 
cus on the role of motility on species interactions within 
a biofilm [50] . 

An aspect that has not been addressed in our work by 
the mean field approach is the nature and role of density 
fluctuations. For the case of a single particle species it has 
been demonstrated that large fluctuations in the density 
can occur [30, 3Tl|33l|5T]. These have been termed giant 
number fluctuations in the context of nematic particle 
interactions Questions that arise are, for example, 
how the nature of these fluctuations changes in the case of 
binary mixtures of self-propelled particles, whether and 
how the two particle species interact through density fluc- 
tuations, and how the situation changes when different 
alignment rules apply. 

Another issue concerns the macroscopic equations de- 
rived in the last section. As discussed, these are stable 
only for small and moderate particle densities. On the 
one hand, it will be interesting to flnd closure relations 
that stabilize these equations also for higher densities, 
even if the deviations from the initial equations increase. 
On the other hand, the non-equilibrium generalization of 
other transformations from particle to fleld descriptions 
that are suitable in the high-density limit is a compelling 
task for the future. 

In conclusion, here, as a flrst step, we have shown nu- 
merical results for the simplest case where the particles 
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of different species feature the same behavior. That is, 
single particles of different species propel with the same 
velocity, show the same orientational diffusion, and fol- 
low the same orientational ordering rules. Rich behavior 
is to be expected when these confining conditions are 
weakened. Investigations for particle species of different 
velocities and different magnitudes of orientational diffu- 
sion are currently underway. 
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